home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / pfe.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  4.2 KB  |  121 lines  |  [TEXT/????]
  1. ;;; $Header: pfe.scm,v 1.5 88/07/18 22:46:59 GMT hal Exp $
  2. ;;;; Partial-fractions expansions
  3.  
  4. (if-mit
  5.  (declare (usual-integrations = + - * /
  6.                  zero? 1+ -1+
  7.                  ;; truncate round floor ceiling
  8.                  sqrt exp log sin cos)))
  9. (define (rat-func->pfe ratfunc)
  10.   (numer-poles->pfe (car ratfunc)
  11.                     (sort (poly->roots (cadr ratfunc))
  12.               (lambda (z1 z2)
  13.                 (< (magnitude z1) (magnitude z2))))))
  14.  
  15. (define (pfe->rat-func pfe)
  16.   (let ((residues (car pfe)) (poles (cadr pfe)))
  17.     (let ((np (pfe->numer-poles residues poles)))
  18.       (list (car np) (roots->poly poles)))))
  19.  
  20.  
  21. ;;; The following partial-fraction procedure returns a list of two terms
  22. ;;; -- (residue-list pole-list) given a list of numer-poly coefficients
  23. ;;; and the denominator pole-list.  For best accuracy, denominator roots
  24. ;;; should be entered in increasing order.  Repeated roots must be
  25. ;;; clustered together.  As an example consider
  26. ;;;
  27. ;;; 2x^5 + 9x^4 - x^3 - 26x^2 + 5x - 1
  28. ;;; ----------------------------------
  29. ;;;     (x + 1)^2 (x - 1)^3 (x + 2)   
  30. ;;;
  31. ;;;                3      -2      -1         1       3     1
  32. ;;;           = ------- + --- + ------- + ------- + --- + ---
  33. ;;;             (x+1)^2   x+1   (x-1)^3   (x-1)^2   x-1   x+2
  34. ;;;
  35. ;;;         (numer-poles->pfe '(-1 5 -26 -1 9 2) '(-1 -1 1 1 1 -2))
  36. ;;;               ==> ((3 -2 -1 1 3 1) (-1 -1 1 1 1 -2))
  37.  
  38. (define (numer-poles->pfe numer-poles)
  39.   (let ((numer-poly (car numer-poles))
  40.     (pole-list (cadr numer-poles)))
  41.     (define (pfe-iter poly poles)
  42.       (if (null? poles)
  43.       '()
  44.       (let ((den (other-roots->poly (cdr poles) (car poles))))
  45.         (let ((coef (/ (poly-value poly (car poles))
  46.                (poly-value den (car poles)))))
  47.           (cons coef
  48.             (pfe-iter
  49.              (car
  50.               (div-polys
  51.                (sub-polys poly (scale-poly coef den))
  52.                (list (- (car poles)) 1)))
  53.              (cdr poles)))))))
  54.     (define (other-roots->poly root-list root)
  55.       (if (member root root-list)
  56.       (other-roots->poly (cdr root-list) root)
  57.       (roots->poly root-list)))
  58.     (list (pfe-iter numer-poly pole-list) pole-list)))
  59.  
  60. ;;; The following procedure reconstructs a rat-func in the form
  61. ;;; (numer-poly pole-list) from a partial-fraction expansion in the
  62. ;;; form (residue-list pole-list). Example
  63. ;;;
  64. ;;;       (pfe->numer-poles '((3 -2 -1 1 3 1) (-1 -1 1 1 1 -2)))
  65. ;;;           ==> ((-1 5 -26 -1 9 2) (-1 -1 1 1 1 -2))
  66.  
  67. (define (pfe->numer-poles pfe)
  68.   (let ((residue-list (car pfe)) (pole-list (cadr pfe)))
  69.     (define (iter residues poles pending-poles)
  70.       (if (null? residues)
  71.       '()
  72.       (add-polys (scale-poly (car residues) (roots->poly (cdr poles)))
  73.              (if (member (car poles) (cdr poles))
  74.              (iter (cdr residues)
  75.                    (cdr poles)
  76.                    (cons (car poles) pending-poles))
  77.              (iter (cdr residues)
  78.                    (append (cdr poles)
  79.                        (cons (car poles) pending-poles))
  80.                    '())))))
  81.     (list (iter (reverse residue-list) (reverse pole-list) '()) pole-list)))
  82.  
  83. ;;; The following procedure returns a function of a real variable, t,
  84. ;;; whose value is the Inverse Laplace Transform of a rational function
  85. ;;; described by a RESIDUE-LIST of residues and a POLE-LIST of poles, as
  86. ;;; might be returned by RAT-FUNC->PFE above.
  87.  
  88. (define (pfe->time pfe)
  89.   (define (rptd-pole-residues poly residue-list pole-list)
  90.     ((if (or (null? (cdr residue-list))
  91.              (not (= (car pole-list) (cadr pole-list))))
  92.      list
  93.      rptd-pole-residues)
  94.      (cons (car residue-list) poly)
  95.      (cdr residue-list)
  96.      (cdr pole-list)))
  97.   (define (insert-factorials lst n)
  98.     (if (null? lst)
  99.     '()
  100.     (cons (car lst)
  101.           (map (lambda (x) (/ x n))
  102.            (insert-factorials (cdr lst) (1+ n))))))
  103.   (let loop ((residue-list (car pfe)) (pole-list (cadr pfe)))
  104.     (if (null? residue-list)
  105.         0
  106.         (let ((proc
  107.            (lambda (t)
  108.          (if (negative? t)
  109.              0
  110.              (exp (* (car pole-list) t))))))
  111.       (if (or (null? (cdr residue-list))
  112.           (not (= (car pole-list) (cadr pole-list))))
  113.           (+ (* (car residue-list) proc)
  114.              (loop (cdr residue-list) (cdr pole-list)))
  115.           (let ((new-lists (rptd-pole-residues (list (car residue-list))
  116.                            (cdr residue-list)
  117.                            (cdr pole-list))))
  118.             (+ (* (poly->procedure (insert-factorials (car new-lists) 1))
  119.               proc)
  120.            (loop (cadr new-lists) (caddr new-lists)))))))))
  121.